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Towards a determination ofc$w using NSPT 



C. Torrero 



1. Motivation 

A problem one has often to face when handling lattice results is taking the continuum limit. 

A viable way of reducing the impact of lattice artifacts (by removing some of them) is given 
by the Symanzik improvement programme [jl|] which, as is well-known, has allowed to remove 
6 [a) artifacts in unquenched simulations. The irrelevant term to be added to the lattice action was 
determined by Sheikoleslami and Wohlert [||] and contains the so-called csw coefficient which can 
be expanded perturbatively in even powers of the bare coupling go- 
While the zero- and one-loop coefficients of this expansion have already been computed for diffe- 
rent lattice actions [^, |]] (see also references therein), the two-loop contribution is still unknown: 
within the Wilson formulation of Lattice QCD (LQCD), we address its determination using NSPT, 
a tool that allows for perturbative calculation in lattice-regularized quantum field theories. 

In the second part, we discuss a technical issue of NSPT that is related to the discretization of 
the Langevin equation that governs the evolution of the system: our target is to obtain high-order 
integrators, aiming at reducing the computer-time at fixed accuracy. 



2. Wilson formulation of LQCD and improvement: some notation 

The Wilson action Sw of LQCD can be decomposed into gauge (Sq) and fermionic (Sp) parts 



with 1 



Sg = /3£ (l-£[Up V (n) + Ut v (n)]), 

fi>v 



(2.1) 



where /3 = 2N c /gQ, N c is the number of colours, U^ v (n) is the lattice plaquette and, 

S F = Vab{n)JZ n ab,mF>c[UWf>c{m) , 

n, a,b 
m,j},c 



(2.2) 



with 



ab.mfic 



4e 



be °~n+p.,m be $*-A, 



+ 



+ (M + 4r)8 nm 8 a p8 l 



'be ■ 



(2.3) 



In eq. (2.3), r is the Wilson parameter (which we will set to 1) while Mq is the bare mass. 

The Sheikoleslami-Wohlert irrelevant contribution (Ssw) to be added to the Wilson action is 
given by 2 , 

Ssw = j c swY, V / («) (J mv^vWV / "(«) > 

n, H,v 

where a MV = i/2[Yn,Yv] while F^ v (n) reads, 



(2.4) 



'in the equations of this and the subsequent section it is understood that all dimensionful quantities have been 
rescaled with powers of the lattice spacing a to be dimensionless (some of them carry an extra "hat" to emphasize this). 

2 In the following equations, spin and colour subscripts are suppressed to ease the notation: they can obviously be 
restored as in eq. (2.3). 
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Fnv(n) = -(Q flv (n)-Q Vfl (n)), (2.5) 

with Qn V (n) being the clover term, Q^ v {n) = U^ v (n) + U^ v ,^(n) +£/ V) _j 1 (n) + U- ll - v {n). 
As already anticipated, the coefficient c$w appearing in eq. (2.4) can be decomposed as 



- _J°) , ilj , J 2 ) „4 



6 ) 



(2.6) 



where cLi is the target of our computation. 



3. How to measure c 



(2) 

sw 



,(2) 



A suitable observable to determine cLp is given by the pion propagator 

Gb c (n — m 

a, p, S, e 



£ (¥ab(n)(r 5 ) a pVl3b{fi)\fr Sc (m)(Y5)se¥ec(m) 

a, p, S, e 



L ((75)«(3 
a, P, S, e 



nflb, mSc 



M- 1 



mec, nab 



a, e 



M 



mec, nab 



M 



mec, nab 



a, e 



Af" 



mec, nab 



(3.1) 



where M is the fermionic operator obtained by adding together eqs. (2.2) and (2.4): for details 
about its inversion, see After switching to momentum space and defining the dimensionless 
quantities = p^a and p 2 = E/zPu (being the p^'s the lattice momentum components), one can 
invert the propagator to obtain the f -function which can be decomposed as, 



r(p,m cng0 ) =p 2 +Mt+M w (p)-Z(p,mcngo) , 



2/A\ 



(3.2) 



where Mw(p) is the irrelevant Wilson mass, Mq the (perturbative) pion rest mass and t,(p,m cr ,go) 
the self-energy with m cr the critical mass defined as m cr = L(0,m cr ,go). Since we want to develop 
a mass-independent scheme, we will both set Mq equal to zero and subtract the proper mass coun- 
terterms to keep fermions massless. 

Given that we are eventually interested in a perturbative approach, we can expand L(p,m cr ,go) 
as a series in gj h i.e. £(/>,ra cr ,go) = Hk^^iP^crjg 1 ^, and decompose a generic coefficient 

(p,m cr ) by means of hypercubic invariants as, 



father) = of (m cr ) + of (m cr ) £pj + af (m cr )£pj + . . . . 



,(*) 



(3.3) 



(2) 

A possible approach to determine c Jn would consist of expanding the pion and quark self 



), then take this value 



energies 4 in a combined way: one could take cf^ = 1 to obtain K c to ^(j3 _1 
to tune Cyjy to make the pion mass vanish; next, one revisits the quark propagator to determine K c 
to ^(j8~ 2 ) and so on till Cy™ is determined 5 . 



3 The subscript "U" means that the corresponding average has to be performed on gauge configurations only. 
4 Formulae similar to eqs. (3.2)-(3.3) hold also for the quark propagator though the Dirac structure is more involved. 
5 It is clear that, requiring the pion to be massless, also implies setting Mq equal to zero. 
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An alternative strategy could be the following: recall that Ssw was introduced to remove 
0(a) artifacts and the pion propagator contains the product of two quark propagators. One should 
be able to establish a correspondance between terms proportional to a in the operator M~ 1 and the 

7 ~ (k) (Q) (I) 

ones proportional to a in the T-function, namely . If one now tunes and c s ^ to their known 

(l) (2) (2) (3) 

values and observes that, correspondingly, ocj and a 3 vanish, one can fix c s ^ by requiring (Xj 
to be zero. This approach is maybe less rigorous but nonetheless should be worth studying. 



4. Numerical setup 

The method of our choice is NSPT. It is related to Stochastic Quantization [^] which consists 
of introducing an extra coordinate, a stochastic time t, together with an evolution equation of the 
Langevin type, 

d(j>(x,t) dS[(j)} 

where in this example <p(x,t) is a scalar field while t]{x,t) is a Gaussian noise. 

Starting from this, the usual Feynman-Gibbs integration can be recovered by noise-averaging as 

Z- l j [D<j)]0{<l,(x)]e- s ^ =Kmjfdt'(0[<l>r,(x,t')}) 7} . (4.2) 

For SU (3) lattice variables the Langevin equation has to be modified into 

d t U^n,t) = -iT A (v n ^ jA S[U] + ri^(n,t))u^n,t) , (4.3) 



in order to obtain an evolution of the variables within the group: here T A = X A /2 are Gell-Mann 
matrices while ri A (n,t) are again Gaussian noise components. 

The missing ingredient, i.e. Perturbation Theory, is introduced by expanding the £/'s as 6 , 

£/ M M^£jHt/«(x,f), (4.4) 
k 

When plugging this into the Langevin equation, this results in a system of coupled differential 
equations that can be solved numerically via a discretization of the stochastic time t = Nz, where 
T is a time step. 

In practice, the system is evolved for different values of T, then we average over each therma- 
lized signal to realize the limit t — > oo of eq. (4.2). Finally we extrapolate in z to the z = limit 
of the desired observable: this extrapolation is needed since the correct Boltzmann equilibrium 
distribution is recovered only for continuous t. 



5. High-order integrators for NSPT 

The csiy-simulations are an ongoing project such that the numerical results presented here 
refer to improvements of the NSPT algorithm. 

6 The expansion in the computer code is thought on jS -1 / 2 rather than go: converting the corresponding coefficients 
is obviously straightforward. 
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As mentioned above, one has to perform simulations with different values of z to extrapolate 
towards the limit z — > and this increases the required computer-time. Since the smaller z is, the 
more iterations N are needed, a possible way to save computer-time might consist of employing 
larger values of z; however, this would compromise the accuracy of the subsequent T-extrapolation. 
The solution to this drawback is well-known and is represented by high-order integrators for the 
Langevin equation: these indeed increase the power of the leading T-dependence thus allowing to 
safely recover the z — > limit even at large values of the time step. 

The easiest way of determining high-order integrators is probably by generalizing the usual 
Runge-Kutta schemes for the scalar case: there, given a scalar variable y(z), its derivative / = 
f(z,y) and an initial value y(zo) = yo, the m-th order integrator reads: 



(5.1) 



y n +\ =y n + ^Y^ blkl {ki=fhn + ciZ,y n +zY^ai, r k r j ; h=f(z n ,y n 

1=1 V V r =l ' 

The generalization to non-Abelian variables appears straightforward: 

m r m . . -, 

yn+l = y n + ^Y^ blkl > ^( X ' T «+l)= ex P -^^/( 7 7m( X ' T «) + ^) U li( x > T n), (5-2) 

1=1 1 1=1 

/ 1 1 \ 



I- 1 

r=l ' A 

where 5 [I/O] is the expression of the action where all variables have been replaced as 

r=l 



(5-3) 



U^(x,z n ) — >exp -(T^a/J^^Tnj+U U^(x,z n 



(5.4) 



where it is understood that ki = £ A T A V Xtjl A S[U (z n )]. 

As is manifest, the number of operations per update increases with the order of the integrator: one 
will eventually be able to employ larger time steps in the simulations, thus reducing the number of 
iterations, but at the price of more costly iterations. Our study seems to indicate that overall savings 
of up to a factor of two can still be achieved. 

Now everything boils down at getting the coefficients ai r ,bi,ci in eq. (5.1) but they can easily 
be found in the literature. 

At present, the highest integrator available for the NSPT Langevin equation is second-order [[7|] 
and reads, 



C/ m (jc,t„ + i) = exp 



'11 + ^jf) Q^i + \*h ) - ' v' *rin (-v. t„ ) 



U[i(x,z n ) 



£r A v,, M , A 5[t/(T n )] , 

A 

Z,T A V x ^ A S[U i2) ] , 

A 

exp - izki - iVzrjn (x, z n ) f/ M (x, z n ) 



(5.5) 
(5-6) 
(5.7) 
(5.8) 
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where Ca is the Casimir invariant of the Lie group's adjoint representation: note that the noise 
employed in eqs. (5.5) and (5.8) is the same. The second term in the first square brackets in 
eq. (5.5) — not appearing in Runge-Kutta literature — comes from the non-commutativity of group 
derivatives introduced in eq. (4.3). Apart from a rescaling 7 of x with /3 [§], the integrator can be 
determined by calculating the equilibrium distribution of the corresponding discretized Fokker- 
Planck equation. 

At present, analytical calculations are undertaken in order to compute the corresponding non- 
Abelian shifts for both the third- and the fourth-order integrator. 

6. Preliminary results 

As can be seen from eq. (5.5), the above-mentioned non-Abelian shift only affects loops higher 
than the first and thus we are already in the position of comparing one-loop results obtained with 
different integrators: the variable we choose is the plaquette at a lattice extent L = 4. 
Figure 1 shows how the slope in x changes with the integrator as expected while, in Table 1, we 
collect our numerical results: as desired, the accuracy remains rather good even when employing 
larger x values in the simulations. 



Order of integrator 


Employed time steps 


1-loop plaquette 


1 


10, 15, 20 


1.9930(7) 


2 


50, 60, 70 


1.9922(6) 


3 


90, 100, 110 


1.9918(10) 


4 


110, 122, 130 


1.9914(10) 



Table 1: Comparison between 1-loop results for different integrators at L=4. the diagrammatic L = 4 value 
reads 1.9921875.... 

The second-order integrator should already be working to any perturbative order so that, as 
a second test, we can check whether higher loops are under control too when increasing the time 
steps in the simulations: this is done in Table 2 where benchmark plaquette results at L = 4 are 
provided by the first-order integrator that has been in use since long. 



Order of integrator 


1st loop 


2nd loop 


3rd loop 


4th loop 


1 


1.9930(7) 


1.2027(18) 


2.8781(67) 


8.994(30) 


2 


1.9922(6) 


1.2002(17) 


2.8778(62) 


8.990(28) 



Table 2: T — * results from 1st and 2nd order integrators at L = 4: the 1- and 2-loop diagrammatic values 
read 1.9921875 . . . and 1.2037037 . . .. The time steps employed can be found in Table 1. 

In Figure 2 we compare the x — >0 results from the second order integrator to the corresponding 
diagrammatic results for different L. At L = 2 there is some disagreement which might be due to 
different ways of treating zero modes. The ratios of diagrammatic finite over infinite volume results 

7 The j8 prefactor within eq. (2. 1) needs to be compensated for and, consequently, the rescaling x v- > t//3 is required 
in every integration scheme, including Euler's most trivial one. 
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Figure 1: 1-loop plaquette vs. T at L=4: data 
come from first-, second- and third-order integra- 
tor (blue, red and black points respectively). 
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Figure 2: Ratio NSPT results/diagrammatic values 
vs. L for the 1- and 2-loop plaquette (blue dots and 
red diamonds respectively). 



read 0.907 (0.907) for L = 2 and 0.994 (0.986) for L = 4 at 1-loop (2-loop) level, respectively: finite 
volume effects are much bigger than this disagreement between diagrammatic results (neglecting 
zero modes) and NSPT (subtracting zero modes). The infinite volume limit remains unaffected. 



7. Conclusions 



(2) 

The computation of c^, is at an early stage: at present, we are trying to single out the most 
reliable approach. 

As for higher-order integrators for the Langevin equation, first results seem to confirm good 
gains in computer-time, without any loss in numerical accuracy. 
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